fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

ex12

hierarchical model  
class c=1~k 
b0c~N(b00,sb0)  
b1c~N(b10,sb1)  
yc~N(b0c+b1c*x,s)
n=50
k=9
x=runif(n,0,20)
c=sample(letters[1:k],n,replace=T)
b00=rnorm(k,10,5)
b0=b00[factor(c)]
b10=rnorm(k,2,1)
b1=b10[factor(c)]
y=rnorm(n,b0+b1*x,5)
qplot(x,y,shape=c,size=I(2))+
  scale_shape_manual(values=1:k)

qplot(x,y,col=c)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
estimate as no class
y~N(b00+b10*x,s)

ex8-0.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
mdl=cmdstan_model('./ex8-0.stan') 
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
fn(mdl,data)
## Initial log joint probability = -7288.33 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       20      -135.446   0.000294158     0.0039898           1           1       35    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__    -135.45
##    b0         7.07
##    b1         1.92
##    s          9.11
##    m0[1]     14.63
##    m0[2]     15.39
##    m0[3]     39.30
##    m0[4]     16.67
##    m0[5]     44.47
##    m0[6]     44.65
##    m0[7]     29.04
##    m0[8]     13.85
##    m0[9]     38.77
##    m0[10]    29.90
##    m0[11]     8.42
##    m0[12]    44.15
##    m0[13]    44.32
##    m0[14]    33.87
##    m0[15]    24.36
##    m0[16]    26.71
##    m0[17]    17.59
##    m0[18]    20.37
##    m0[19]    10.98
##    m0[20]    15.36
##    m0[21]    11.53
##    m0[22]    36.15
##    m0[23]    40.70
##    m0[24]    15.24
##    m0[25]    23.22
##    m0[26]    18.09
##    m0[27]    36.32
##    m0[28]    23.05
##    m0[29]    30.48
##    m0[30]    36.81
##    m0[31]    19.11
##    m0[32]    22.81
##    m0[33]    17.16
##    m0[34]    43.25
##    m0[35]    35.80
##    m0[36]    19.27
##    m0[37]    19.67
##    m0[38]    37.20
##    m0[39]    32.30
##    m0[40]    34.85
##    m0[41]    31.38
##    m0[42]    41.09
##    m0[43]    13.92
##    m0[44]    27.25
##    m0[45]    40.50
##    m0[46]    41.42
##    m0[47]    25.03
##    m0[48]    20.34
##    m0[49]    19.14
##    m0[50]    41.19
##    y0[1]     12.35
##    y0[2]      7.84
##    y0[3]     40.86
##    y0[4]     19.76
##    y0[5]     57.44
##    y0[6]     48.30
##    y0[7]     39.57
##    y0[8]     14.43
##    y0[9]     30.93
##    y0[10]    28.43
##    y0[11]    -5.66
##    y0[12]    41.27
##    y0[13]    43.45
##    y0[14]    34.74
##    y0[15]    23.13
##    y0[16]     9.36
##    y0[17]    17.71
##    y0[18]    30.73
##    y0[19]     5.88
##    y0[20]     7.10
##    y0[21]    -0.90
##    y0[22]    38.26
##    y0[23]    44.15
##    y0[24]     3.58
##    y0[25]    13.92
##    y0[26]    26.23
##    y0[27]    38.65
##    y0[28]    10.51
##    y0[29]    22.36
##    y0[30]    34.38
##    y0[31]    45.13
##    y0[32]    28.55
##    y0[33]     8.98
##    y0[34]    54.68
##    y0[35]    39.62
##    y0[36]    27.56
##    y0[37]    18.59
##    y0[38]    31.25
##    y0[39]    24.49
##    y0[40]    37.17
##    y0[41]    35.15
##    y0[42]    47.67
##    y0[43]    23.43
##    y0[44]    25.27
##    y0[45]    25.17
##    y0[46]    46.84
##    y0[47]    17.49
##    y0[48]    31.01
##    y0[49]    22.22
##    y0[50]    36.26
##    m1[1]      8.42
##    m1[2]     12.45
##    m1[3]     16.47
##    m1[4]     20.50
##    m1[5]     24.52
##    m1[6]     28.55
##    m1[7]     32.57
##    m1[8]     36.60
##    m1[9]     40.62
##    m1[10]    44.65
##    y1[1]      1.69
##    y1[2]     -3.07
##    y1[3]     13.57
##    y1[4]     17.58
##    y1[5]     27.78
##    y1[6]     31.94
##    y1[7]     20.63
##    y1[8]     48.37
##    y1[9]     26.64
##    y1[10]    50.46
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable    mean  median    sd   mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   -134.70 -134.43  1.18  0.96 -136.92 -133.39 1.00      654     1121
##    b0        7.27    7.27  2.65  2.73    2.96   11.49 1.01      353      302
##    b1        1.90    1.90  0.22  0.23    1.52    2.25 1.00      435      606
##    s         9.54    9.43  1.02  0.99    8.04   11.30 1.00     1851     1491
##    m0[1]    14.76   14.77  1.94  2.01   11.61   17.93 1.01      372      339
##    m0[2]    15.51   15.56  1.88  1.94   12.45   18.58 1.01      377      344
##    m0[3]    39.19   39.22  1.98  1.99   35.88   42.33 1.00     1414     1249
##    m0[4]    16.77   16.80  1.77  1.81   13.86   19.69 1.01      389      344
##    m0[5]    44.31   44.31  2.46  2.47   40.30   48.34 1.00      994     1094
##    m0[6]    44.49   44.49  2.48  2.50   40.46   48.56 1.00      984     1106
##    m0[7]    29.04   29.06  1.35  1.34   26.79   31.24 1.00     1560     1204
##    m0[8]    13.99   14.00  2.01  2.07   10.74   17.23 1.01      367      337
##    m0[9]    38.67   38.70  1.93  1.94   35.48   41.74 1.00     1474     1259
##    m0[10]   29.88   29.90  1.38  1.37   27.60   32.10 1.00     1763     1238
##    m0[11]    8.61    8.60  2.52  2.57    4.55   12.62 1.01      354      296
##    m0[12]   44.00   44.01  2.43  2.44   40.03   47.98 1.00     1011     1083
##    m0[13]   44.17   44.17  2.45  2.46   40.18   48.18 1.00     1001     1094
##    m0[14]   33.81   33.85  1.57  1.61   31.25   36.40 1.00     2017     1399
##    m0[15]   24.40   24.39  1.36  1.31   22.18   26.68 1.00      684      815
##    m0[16]   26.72   26.73  1.33  1.30   24.53   28.94 1.00     1011      961
##    m0[17]   17.69   17.71  1.70  1.73   14.90   20.48 1.01      400      358
##    m0[18]   20.44   20.45  1.52  1.51   17.87   22.97 1.00      457      447
##    m0[19]   11.14   11.15  2.27  2.31    7.49   14.75 1.01      357      311
##    m0[20]   15.48   15.53  1.88  1.94   12.42   18.55 1.01      377      349
##    m0[21]   11.68   11.69  2.22  2.27    8.12   15.21 1.01      358      305
##    m0[22]   36.07   36.13  1.72  1.76   33.25   38.83 1.00     1752     1376
##    m0[23]   40.59   40.59  2.11  2.12   37.11   43.97 1.00     1260     1177
##    m0[24]   15.36   15.40  1.89  1.95   12.30   18.45 1.01      376      344
##    m0[25]   23.27   23.28  1.39  1.36   20.95   25.64 1.00      589      720
##    m0[26]   18.19   18.20  1.67  1.67   15.43   20.95 1.01      407      396
##    m0[27]   36.25   36.31  1.74  1.77   33.41   39.01 1.00     1734     1276
##    m0[28]   23.10   23.11  1.40  1.38   20.76   25.48 1.00      578      719
##    m0[29]   30.45   30.46  1.40  1.41   28.15   32.72 1.00     1793     1254
##    m0[30]   36.73   36.80  1.78  1.79   33.84   39.56 1.00     1677     1319
##    m0[31]   19.20   19.19  1.60  1.60   16.54   21.80 1.01      426      412
##    m0[32]   22.86   22.87  1.40  1.39   20.50   25.26 1.00      563      703
##    m0[33]   17.27   17.28  1.74  1.78   14.43   20.11 1.01      395      358
##    m0[34]   43.10   43.11  2.34  2.37   39.24   46.94 1.00     1065     1082
##    m0[35]   35.73   35.78  1.70  1.74   32.95   38.45 1.00     1791     1354
##    m0[36]   19.35   19.35  1.59  1.58   16.71   21.95 1.01      429      423
##    m0[37]   19.75   19.76  1.56  1.56   17.13   22.34 1.01      439      432
##    m0[38]   37.12   37.17  1.81  1.81   34.16   40.00 1.00     1678     1302
##    m0[39]   32.26   32.28  1.48  1.52   29.80   34.67 1.00     2022     1413
##    m0[40]   34.79   34.83  1.63  1.69   32.14   37.42 1.00     1888     1422
##    m0[41]   31.35   31.35  1.43  1.46   28.99   33.69 1.00     1910     1288
##    m0[42]   40.97   40.98  2.14  2.15   37.45   44.41 1.00     1226     1177
##    m0[43]   14.05   14.07  2.00  2.07   10.82   17.29 1.01      367      337
##    m0[44]   27.25   27.26  1.33  1.30   25.07   29.46 1.00     1121     1083
##    m0[45]   40.39   40.40  2.09  2.11   36.94   43.74 1.00     1294     1177
##    m0[46]   41.29   41.30  2.17  2.19   37.72   44.78 1.00     1199     1177
##    m0[47]   25.06   25.05  1.34  1.31   22.85   27.33 1.00      755      875
##    m0[48]   20.42   20.43  1.52  1.51   17.84   22.95 1.00      456      447
##    m0[49]   19.23   19.23  1.60  1.59   16.58   21.83 1.01      426      412
##    m0[50]   41.07   41.08  2.15  2.16   37.53   44.52 1.00     1217     1177
##    y0[1]    14.68   15.02  9.79  9.53   -1.73   30.53 1.00     1621     1929
##    y0[2]    15.05   15.16  9.65  9.51   -0.69   30.65 1.00     1965     1745
##    y0[3]    39.40   39.43  9.59  9.26   23.70   55.27 1.00     1975     1636
##    y0[4]    16.81   16.82  9.87  9.70    0.76   32.58 1.00     1823     1853
##    y0[5]    44.67   45.01 10.13 10.08   27.97   60.55 1.00     1899     2013
##    y0[6]    44.64   44.56  9.88  9.68   28.74   60.54 1.00     1778     1625
##    y0[7]    29.20   29.06  9.56  9.33   13.46   45.18 1.00     2004     1993
##    y0[8]    14.53   14.83  9.81  9.82   -1.54   31.22 1.00     1718     1730
##    y0[9]    38.93   39.14  9.80  9.91   23.28   55.07 1.00     2074     2147
##    y0[10]   30.01   30.33  9.54  9.16   14.32   45.85 1.00     2121     1897
##    y0[11]    8.88    8.83 10.01 10.16   -7.32   25.09 1.00     1425     1765
##    y0[12]   44.11   43.86  9.78  9.63   28.26   60.57 1.00     1703     1968
##    y0[13]   44.60   44.48 10.03 10.00   28.72   61.03 1.00     1907     1914
##    y0[14]   33.62   33.67  9.46  9.52   18.11   48.91 1.00     2084     1776
##    y0[15]   24.32   24.40  9.75  9.29    8.05   40.21 1.00     2011     1927
##    y0[16]   26.52   26.66  9.89  9.71   10.17   42.52 1.00     1985     1895
##    y0[17]   17.40   17.49  9.40  9.11    2.24   33.05 1.00     1467     1856
##    y0[18]   20.17   20.37  9.76  9.50    3.65   35.48 1.00     1773     1918
##    y0[19]   10.92   11.05 10.03  9.74   -6.40   27.30 1.00     1605     1822
##    y0[20]   15.22   15.09  9.62  9.28   -0.21   31.71 1.00     2012     2040
##    y0[21]   11.84   11.75  9.78  9.40   -4.31   27.91 1.00     2007     2022
##    y0[22]   35.99   36.04  9.85  9.54   19.75   52.74 1.00     1839     1966
##    y0[23]   40.30   40.20  9.73  9.50   24.32   56.55 1.00     1850     2101
##    y0[24]   14.98   15.16  9.78  9.72   -1.03   30.74 1.00     1740     1933
##    y0[25]   23.10   23.08  9.69  9.65    7.18   38.83 1.00     1985     1894
##    y0[26]   17.90   17.87  9.56  9.44    2.39   33.49 1.00     1725     1820
##    y0[27]   36.20   35.89 10.09  9.99   20.17   52.61 1.00     1887     1857
##    y0[28]   23.20   23.17  9.79  9.75    7.14   39.26 1.00     2120     1862
##    y0[29]   30.48   30.58  9.69  9.59   14.40   46.51 1.00     1978     2057
##    y0[30]   36.49   36.71  9.74  9.52   19.77   52.15 1.00     1882     1834
##    y0[31]   19.30   19.20  9.85  9.73    3.22   35.19 1.00     1874     2056
##    y0[32]   22.49   22.65  9.68  9.29    6.64   38.45 1.00     1540     1815
##    y0[33]   17.16   17.12  9.71  9.74    1.61   33.53 1.00     1487     1877
##    y0[34]   43.20   42.82  9.93  9.95   27.64   59.43 1.00     1841     2053
##    y0[35]   35.64   36.00  9.82 10.04   19.41   51.69 1.00     1877     1928
##    y0[36]   18.87   18.75  9.46  9.19    3.48   34.64 1.00     1884     2007
##    y0[37]   19.68   20.00  9.65  9.67    3.59   35.53 1.00     2008     1930
##    y0[38]   37.22   37.16  9.75  9.58   21.44   53.11 1.00     2254     1746
##    y0[39]   32.45   32.55  9.47  9.06   16.60   47.96 1.00     2158     2028
##    y0[40]   34.78   34.59 10.04  9.82   18.71   51.53 1.00     2179     1913
##    y0[41]   31.37   31.31  9.79  9.87   15.72   47.26 1.00     2046     1826
##    y0[42]   41.12   41.15  9.96  9.85   24.75   57.09 1.00     1966     1957
##    y0[43]   13.91   14.12  9.75  9.50   -1.81   29.72 1.00     1893     1929
##    y0[44]   27.41   27.43  9.72  9.68   11.01   43.45 1.00     2043     1957
##    y0[45]   40.08   39.67  9.84  9.63   24.56   56.34 1.00     1896     2101
##    y0[46]   40.98   40.70  9.51  9.17   25.41   57.13 1.01     1803     2056
##    y0[47]   25.19   25.18  9.68  9.65    9.21   40.73 1.00     1930     1845
##    y0[48]   20.12   20.13  9.69  9.52    4.25   36.20 1.00     1914     2003
##    y0[49]   19.14   19.28  9.71  9.85    2.66   35.18 1.00     1848     2053
##    y0[50]   41.18   41.06  9.84  9.92   25.22   57.56 1.00     1907     1636
##    m1[1]     8.61    8.60  2.52  2.57    4.55   12.62 1.01      354      296
##    m1[2]    12.60   12.62  2.13  2.19    9.15   16.04 1.01      361      310
##    m1[3]    16.58   16.62  1.79  1.82   13.64   19.53 1.01      387      344
##    m1[4]    20.57   20.58  1.51  1.50   18.01   23.09 1.00      461      450
##    m1[5]    24.56   24.55  1.35  1.32   22.34   26.83 1.00      700      802
##    m1[6]    28.55   28.56  1.34  1.33   26.32   30.75 1.00     1438     1209
##    m1[7]    32.53   32.53  1.49  1.53   30.05   34.96 1.00     2025     1422
##    m1[8]    36.52   36.60  1.76  1.79   33.65   39.32 1.00     1701     1276
##    m1[9]    40.51   40.51  2.10  2.12   37.04   43.88 1.00     1283     1177
##    m1[10]   44.49   44.49  2.48  2.50   40.46   48.56 1.00      984     1106
##    y1[1]     8.90    8.97 10.10  9.99   -7.86   25.42 1.00     1445     1904
##    y1[2]    12.67   12.83  9.86  9.54   -3.12   28.87 1.00     1660     1618
##    y1[3]    16.78   16.90  9.85  9.65    0.21   32.38 1.00     1888     1563
##    y1[4]    20.53   20.20  9.79  9.65    4.45   36.95 1.00     1799     1936
##    y1[5]    24.22   24.12  9.68  9.40    8.24   40.32 1.00     1875     1854
##    y1[6]    28.47   28.22  9.73  9.96   12.88   44.60 1.00     2045     1845
##    y1[7]    32.72   32.69  9.54  9.57   16.36   48.27 1.00     1859     1974
##    y1[8]    36.74   36.83  9.57  9.49   21.17   52.02 1.00     1462     2097
##    y1[9]    40.22   40.23 10.00 10.23   23.96   56.12 1.00     1823     1750
##    y1[10]   44.40   44.12  9.71  9.45   28.55   60.17 1.00     1662     1933

estimate as independent class  
all b0l,b1l do not have a distribution  
class c=1~k 
yc~N(b0c+b1c*x,s)  

ex12-1.stan

data{
  int N;
  int K;
  vector[N] x;
  vector[N] y;
  array[N] int c;
}
parameters{
  vector[K] b0;
  vector[K] b1;
  real<lower=0> s;
}
model{
  for(i in 1:N)
    y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0[c[i]]+b1[c[i]]*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
}
mdl=cmdstan_model('./ex12-1.stan') 
data=list(N=n,K=k,x=x,y=y,c=factor(c))

mle=mdl$optimize(data=data)
## Initial log joint probability = -1.23086e+06 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmpwb3xsj/model-20a93d7d50fa.stan', line 15, column 4 to column 42) 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmpwb3xsj/model-20a93d7d50fa.stan', line 15, column 4 to column 42) 
## Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/Rtmpwb3xsj/model-20a93d7d50fa.stan', line 15, column 4 to column 42) 
## Error evaluating model log probability: Non-finite gradient. 
##       99      -87.0685      0.426789       1.49779           1           1      166    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      199      -86.9818    0.00250448     0.0405116      0.6326      0.6326      275    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      222      -86.9817   0.000339057    0.00138477           1           1      300    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -86.98
##    b0[1]     -1.59
##    b0[2]     15.59
##    b0[3]     21.29
##    b0[4]      5.07
##    b0[5]      5.36
##    b0[6]      3.55
##    b0[7]      8.14
##    b0[8]     -2.70
##    b0[9]     -1.45
##    b1[1]      3.85
##    b1[2]      1.03
##    b1[3]      1.70
##    b1[4]      0.18
##    b1[5]      2.14
##    b1[6]      2.27
##    b1[7]      1.64
##    b1[8]      1.55
##    b1[9]      3.17
##    s          3.45
##    m0[1]     12.49
##    m0[2]     14.65
##    m0[3]     32.99
##    m0[4]     29.80
##    m0[5]     47.14
##    m0[6]     35.88
##    m0[7]     15.03
##    m0[8]     27.30
##    m0[9]     49.41
##    m0[10]    44.29
##    m0[11]     5.15
##    m0[12]    47.43
##    m0[13]    46.98
##    m0[14]    18.92
##    m0[15]    24.01
##    m0[16]    27.30
##    m0[17]    15.98
##    m0[18]     8.03
##    m0[19]     8.18
##    m0[20]    20.06
##    m0[21]     5.48
##    m0[22]    20.76
##    m0[23]    24.44
##    m0[24]    15.15
##    m0[25]    35.61
##    m0[26]    16.59
##    m0[27]    47.24
##    m0[28]    22.46
##    m0[29]    37.34
##    m0[30]    38.59
##    m0[31]    18.50
##    m0[32]    35.25
##    m0[33]     5.44
##    m0[34]    53.38
##    m0[35]    31.10
##    m0[36]     6.19
##    m0[37]    22.39
##    m0[38]    39.03
##    m0[39]    33.55
##    m0[40]    31.98
##    m0[41]    32.52
##    m0[42]    51.47
##    m0[43]    11.65
##    m0[44]    27.42
##    m0[45]    33.64
##    m0[46]    37.61
##    m0[47]    37.22
##    m0[48]    33.06
##    m0[49]     7.04
##    m0[50]    43.92
##    y0[1]     15.49
##    y0[2]     15.04
##    y0[3]     36.53
##    y0[4]     36.96
##    y0[5]     49.63
##    y0[6]     33.23
##    y0[7]     14.31
##    y0[8]     24.11
##    y0[9]     48.09
##    y0[10]    46.05
##    y0[11]     4.25
##    y0[12]    50.60
##    y0[13]    46.70
##    y0[14]    19.58
##    y0[15]    22.96
##    y0[16]    30.07
##    y0[17]    16.09
##    y0[18]     8.47
##    y0[19]    11.16
##    y0[20]    16.61
##    y0[21]     4.92
##    y0[22]    26.14
##    y0[23]    30.23
##    y0[24]    10.91
##    y0[25]    38.87
##    y0[26]    12.09
##    y0[27]    47.69
##    y0[28]    20.29
##    y0[29]    44.66
##    y0[30]    41.58
##    y0[31]    16.53
##    y0[32]    36.94
##    y0[33]     4.66
##    y0[34]    57.22
##    y0[35]    27.04
##    y0[36]     5.57
##    y0[37]    19.46
##    y0[38]    39.37
##    y0[39]    38.38
##    y0[40]    31.87
##    y0[41]    32.49
##    y0[42]    51.40
##    y0[43]    16.02
##    y0[44]    28.10
##    y0[45]    30.18
##    y0[46]    45.54
##    y0[47]    35.36
##    y0[48]    40.43
##    y0[49]     3.77
##    y0[50]    44.06
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 2.4 seconds.
## Chain 2 finished in 2.4 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 2.5 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 2.6 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 2.5 seconds.
## Total execution time: 2.7 seconds.
seeMCMC(mcmc,'m')
##  variable    mean median     sd    mad       q5    q95 rhat ess_bulk ess_tail
##    lp__    -97.13 -96.75   3.75   3.65  -103.98 -91.71 1.01      653     1058
##    b0[1]  -333.81   0.10 729.64 372.80 -1710.25 379.46 1.88        5       28
##    b0[2]    15.71  15.82   4.52   4.37     8.11  23.35 1.00     2422     1268
##    b0[3]    21.20  21.26   3.15   3.08    16.07  26.12 1.00     3419     1349
##    b0[4]     5.07   4.87   7.59   7.32    -7.04  17.87 1.00     1171     1040
##    b0[5]     5.42   5.59   4.90   4.84    -2.85  13.28 1.01     2325     1464
##    b0[6]     3.58   3.59   2.37   2.33    -0.11   7.53 1.01     2681      923
##    b0[7]     8.17   8.08   6.32   6.28    -2.23  18.96 1.01     1957     1206
##    b0[8]    -2.66  -2.68   4.49   4.40    -9.93   4.68 1.00     2217     1172
##    b0[9]    -1.40  -1.74   7.34   7.42   -13.18  10.80 1.00     1379     1164
##    b1[1]    31.73   3.79  61.22  31.14   -28.31 147.29 1.88        5       27
##    b1[2]     1.03   1.03   0.31   0.30     0.50   1.55 1.00     2208     1466
##    b1[3]     1.71   1.70   0.26   0.27     1.29   2.12 1.00     3062     1560
##    b1[4]     0.14   0.18   1.57   1.48    -2.53   2.56 1.01     1208     1130
##    b1[5]     2.14   2.13   0.34   0.33     1.58   2.70 1.00     2358     1338
##    b1[6]     2.27   2.27   0.23   0.23     1.89   2.63 1.00     2213     1566
##    b1[7]     1.63   1.64   0.47   0.46     0.88   2.39 1.01     1750     1379
##    b1[8]     1.54   1.54   0.38   0.37     0.92   2.17 1.00     2291     1259
##    b1[9]     3.17   3.18   0.85   0.85     1.76   4.57 1.00     1548     1251
##    s         4.44   4.38   0.57   0.54     3.62   5.47 1.00     1200     1366
##    m0[1]    12.52  12.53   1.74   1.70     9.67  15.34 1.00     2559     1198
##    m0[2]    14.69  14.82   3.55   3.51     8.74  20.37 1.00     2314     1339
##    m0[3]    33.00  33.00   2.15   2.00    29.53  36.55 1.01     1830      989
##    m0[4]    29.75  29.82   2.08   2.04    26.38  33.04 1.00     3134     1408
##    m0[5]    47.12  47.11   2.53   2.45    43.10  51.29 1.00     2200     1276
##    m0[6]    35.87  35.87   2.71   2.61    31.45  40.30 1.01     1840     1047
##    m0[7]    15.03  15.06   1.69   1.63    12.25  17.84 1.00     2363     1303
##    m0[8]    27.24  27.31   2.37   2.34    23.39  30.93 1.00     3277     1371
##    m0[9]    49.43  49.40   2.07   2.04    46.01  52.82 1.00     2173     1445
##    m0[10]   44.41  44.52   4.43   4.06    37.05  51.77 1.00     2402     1261
##    m0[11]    5.18   5.19   2.24   2.17     1.65   8.86 1.01     2700      942
##    m0[12]   47.45  47.44   3.02   2.98    42.57  52.33 1.00     1986     1567
##    m0[13]   46.96  46.95   2.51   2.43    42.98  51.11 1.00     2198     1276
##    m0[14]   18.91  18.92   2.05   2.00    15.45  22.26 1.00     2295     1404
##    m0[15]   24.03  24.02   1.49   1.47    21.61  26.42 1.00     2085     1807
##    m0[16]   27.32  27.34   1.99   1.96    24.07  30.61 1.00     2224     1675
##    m0[17]   16.03  15.98   3.40   3.52    10.53  21.72 1.00     1359     1331
##    m0[18]    8.04   7.99   2.25   2.20     4.47  11.84 1.00     2303     1460
##    m0[19]    8.21   8.21   2.02   1.96     5.01  11.40 1.00     2716      976
##    m0[20]   20.16  20.28   3.33   3.24    14.64  25.76 1.00     2438     1226
##    m0[21]    5.40   5.40   4.57   4.53    -1.99  13.23 1.00     1246     1275
##    m0[22]   20.75  20.78   2.35   2.30    16.78  24.59 1.00     2296     1266
##    m0[23]   24.42  24.44   3.06   2.95    19.30  29.34 1.00     2327     1141
##    m0[24]   15.14  15.11   4.58   4.55     7.70  22.96 1.00     1990     1351
##    m0[25]   35.59  35.64   1.57   1.49    33.01  38.13 1.00     2578     1589
##    m0[26]   16.62  16.63   1.56   1.50    14.05  19.08 1.00     2337     1647
##    m0[27]   47.26  47.24   1.85   1.80    44.19  50.28 1.00     2115     1419
##    m0[28]   22.49  22.48   1.47   1.42    20.07  24.86 1.00     2144     1874
##    m0[29]   37.40  37.42   4.39   4.30    30.11  44.56 1.00     1883     1386
##    m0[30]   38.58  38.54   1.72   1.60    35.80  41.40 1.00     2142     1226
##    m0[31]   18.56  18.51   3.01   3.12    13.54  23.48 1.00     1418     1346
##    m0[32]   35.22  35.27   1.59   1.51    32.61  37.80 1.00     2612     1598
##    m0[33]    5.46   5.38   2.72   2.66     1.00   9.92 1.00     2257     1370
##    m0[34]   53.42  53.42   2.54   2.51    49.23  57.55 1.00     2292     1555
##    m0[35]   31.13  31.12   1.92   1.82    28.05  34.29 1.01     1937      962
##    m0[36]    5.97   6.12   4.57   4.51    -1.58  13.37 1.01     1444     1220
##    m0[37]   22.47  22.54   2.77   2.69    17.99  27.08 1.00     2497     1237
##    m0[38]   39.02  38.98   1.75   1.64    36.20  41.85 1.00     2140     1211
##    m0[39]   33.55  33.55   1.62   1.54    30.88  36.25 1.00     2219     1535
##    m0[40]   31.88  31.91   2.78   2.78    27.42  36.30 1.01     1935     1388
##    m0[41]   32.53  32.54   1.65   1.56    29.85  35.30 1.00     2210     1587
##    m0[42]   51.50  51.49   2.30   2.28    47.74  55.27 1.00     2239     1525
##    m0[43]   11.68  11.70   1.79   1.75     8.78  14.56 1.00     2622     1164
##    m0[44]   27.45  27.41   1.58   1.64    24.91  30.00 1.00     2042     1560
##    m0[45]   33.65  33.69   2.26   2.13    30.02  37.40 1.01     1822      959
##    m0[46]   37.48  37.47   3.66   3.61    31.28  43.25 1.01     1749     1412
##    m0[47]   37.20  37.24   1.49   1.42    34.75  39.65 1.00     2429     1517
##    m0[48]   33.02  33.07   1.76   1.67    30.15  35.86 1.00     2851     1351
##    m0[49]    7.06   6.97   2.42   2.35     3.15  11.09 1.00     2292     1379
##    m0[50]   43.95  43.94   2.71   2.67    39.60  48.31 1.00     1980     1545
##    y0[1]    12.58  12.59   4.92   4.77     4.42  20.43 1.00     2274     2059
##    y0[2]    14.79  14.78   5.68   5.56     5.09  23.90 1.00     2165     1926
##    y0[3]    32.98  33.05   4.99   4.79    24.71  41.03 1.00     1785     1743
##    y0[4]    29.74  29.71   4.80   4.53    21.71  37.69 1.00     2109     1881
##    y0[5]    47.00  46.98   5.04   5.19    38.68  55.10 1.00     2102     1840
##    y0[6]    35.89  35.78   5.26   5.15    27.38  44.81 1.00     1919     1810
##    y0[7]    15.07  14.95   4.70   4.65     7.39  22.91 1.00     1882     1970
##    y0[8]    27.27  27.38   4.98   4.81    19.17  35.38 1.00     1966     1826
##    y0[9]    49.41  49.32   5.00   4.89    40.94  57.59 1.00     1905     2011
##    y0[10]   44.33  44.47   6.34   6.04    33.97  54.55 1.00     2078     1732
##    y0[11]    5.01   5.05   5.02   4.84    -3.22  13.50 1.00     1903     1669
##    y0[12]   47.33  47.31   5.32   5.28    38.76  55.98 1.00     2159     1931
##    y0[13]   47.02  46.90   5.09   5.08    38.52  55.36 1.00     2087     2039
##    y0[14]   18.82  19.05   4.99   5.03    10.71  26.69 1.00     1812     1883
##    y0[15]   23.83  23.99   4.65   4.46    16.06  31.45 1.00     1620     1820
##    y0[16]   27.36  27.28   4.89   4.69    19.55  35.63 1.00     2189     2013
##    y0[17]   15.93  15.79   5.68   5.83     6.78  25.35 1.00     1751     1656
##    y0[18]    7.99   7.92   5.07   4.86    -0.37  16.54 1.00     2014     1777
##    y0[19]    8.14   7.93   4.95   4.78     0.04  16.28 1.00     2163     1917
##    y0[20]   20.05  20.02   5.59   5.65    10.75  29.32 1.00     2254     2014
##    y0[21]    5.44   5.45   6.50   6.52    -5.18  16.30 1.00     1376     1461
##    y0[22]   20.85  20.90   5.21   5.19    12.21  29.18 1.00     2065     1931
##    y0[23]   24.40  24.47   5.35   5.17    15.66  33.13 1.00     2176     1973
##    y0[24]   15.23  15.29   6.44   6.61     5.02  25.62 1.00     2036     1729
##    y0[25]   35.58  35.47   4.79   4.71    27.63  43.48 1.00     2053     1715
##    y0[26]   16.57  16.49   4.68   4.51     9.08  24.25 1.00     2123     1931
##    y0[27]   47.23  47.24   4.80   4.67    39.50  55.17 1.00     2017     1813
##    y0[28]   22.33  22.35   4.66   4.61    14.45  29.91 1.00     1932     1880
##    y0[29]   37.46  37.37   6.23   6.10    27.43  47.88 1.00     2015     2046
##    y0[30]   38.62  38.75   4.79   4.44    30.58  46.50 1.00     1920     1845
##    y0[31]   18.45  18.44   5.51   5.22     9.52  27.59 1.00     1689     1728
##    y0[32]   35.11  35.01   4.93   4.73    26.92  43.24 1.00     1873     1894
##    y0[33]    5.62   5.80   5.20   5.06    -2.95  14.01 1.00     2202     1836
##    y0[34]   53.39  53.32   5.22   5.10    44.92  62.26 1.00     2037     2056
##    y0[35]   31.20  31.30   4.91   4.90    23.21  39.19 1.00     1814     1932
##    y0[36]    6.11   6.09   6.42   6.24    -3.85  16.83 1.00     1800     1731
##    y0[37]   22.54  22.51   5.18   4.91    13.99  30.85 1.00     1950     1779
##    y0[38]   39.06  38.97   4.76   4.62    31.28  46.82 1.00     1975     2010
##    y0[39]   33.49  33.51   4.76   4.86    25.52  41.37 1.00     2051     1917
##    y0[40]   31.93  32.04   5.33   5.23    23.14  40.63 1.00     2008     1814
##    y0[41]   32.52  32.43   4.72   4.69    25.07  40.61 1.00     2106     1865
##    y0[42]   51.44  51.36   4.82   4.84    43.64  59.17 1.00     1932     1805
##    y0[43]   11.54  11.51   4.79   4.69     3.71  19.50 1.00     2166     1568
##    y0[44]   27.43  27.42   4.80   4.63    19.66  35.34 1.00     1900     1821
##    y0[45]   33.77  33.85   4.98   4.80    25.50  41.94 1.00     2077     1946
##    y0[46]   37.66  37.64   5.77   5.79    28.64  46.88 1.00     1867     1618
##    y0[47]   37.23  37.29   4.69   4.62    29.58  44.94 1.00     1820     2092
##    y0[48]   33.06  32.94   4.94   4.70    25.15  41.17 1.00     2165     1787
##    y0[49]    7.17   7.23   5.05   4.96    -1.29  15.50 1.00     2082     2057
##    y0[50]   43.94  44.00   5.22   5.11    35.67  52.61 1.00     2167     1865

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,shape=c,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  #geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  #geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

qplot(x,y,facets=~c,size=I(2),col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

estimate as class have relation each other  
all b0l,b1l have a distribution  
class c=1~k 
b0c~N(b00,sb0)  
b1c~N(b10,sb1)  
yc~N(b0c+b1c*x,s)

ex12-2.stan

data{
  int N;
  int K;
  vector[N] x;
  vector[N] y;
  array[N] int c;
}
parameters{
  real b00;
  real<lower=0,upper=100> sb0;
  vector[K] b0;
  real b10;
  real<lower=0,upper=100> sb1;
  vector[K] b1;
  real<lower=0,upper=100> s;
}
model{
  b0~normal(b00,sb0);
  b1~normal(b10,sb1);
  for(i in 1:N)
    y[i]~normal(b0[c[i]]+b1[c[i]]*x[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0[c[i]]+b1[c[i]]*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
}
mdl=cmdstan_model('./ex12-2.stan') 
data=list(N=n,K=k,x=x,y=y,c=factor(c))

mle=mdl$optimize(data=data)
## Initial log joint probability = -272.182 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       99      -20.0657      0.276822        604236      0.5672      0.5672      141    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      121      -16.6215   8.53733e-05       290.799           1           1      167    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
try(print(mle))
##  variable estimate
##    lp__     -16.62
##    b00        1.02
##    sb0        0.00
##    b0[1]      1.02
##    b0[2]      1.02
##    b0[3]      1.02
##    b0[4]      1.02
##    b0[5]      1.02
##    b0[6]      1.02
##    b0[7]      1.02
##    b0[8]      1.02
##    b0[9]      1.02
##    b10        2.30
##    sb1        0.20
##    b1[1]      2.72
##    b1[2]      2.05
##    b1[3]      2.96
##    b1[4]      2.03
##    b1[5]      2.36
##    b1[6]      2.61
##    b1[7]      1.52
##    b1[8]      1.37
##    b1[9]      1.93
##    s         13.83
##    m0[1]     11.32
##    m0[2]     11.26
##    m0[3]     35.45
##    m0[4]     15.84
##    m0[5]     47.05
##    m0[6]     41.17
##    m0[7]     16.74
##    m0[8]     11.49
##    m0[9]     49.95
##    m0[10]    33.40
##    m0[11]     2.87
##    m0[12]    51.54
##    m0[13]    46.88
##    m0[14]    20.19
##    m0[15]    24.58
##    m0[16]    25.20
##    m0[17]    11.60
##    m0[18]    10.53
##    m0[19]     6.35
##    m0[20]     9.88
##    m0[21]     5.74
##    m0[22]    21.82
##    m0[23]    25.08
##    m0[24]     7.50
##    m0[25]    25.95
##    m0[26]    16.04
##    m0[27]    46.18
##    m0[28]    22.79
##    m0[29]    24.57
##    m0[30]    37.63
##    m0[31]    13.13
##    m0[32]    25.32
##    m0[33]     8.24
##    m0[34]    56.87
##    m0[35]    31.72
##    m0[36]    13.93
##    m0[37]    14.49
##    m0[38]    38.11
##    m0[39]    32.08
##    m0[40]    23.04
##    m0[41]    30.95
##    m0[42]    53.55
##    m0[43]    10.35
##    m0[44]    28.51
##    m0[45]    36.75
##    m0[46]    28.25
##    m0[47]    28.75
##    m0[48]    21.51
##    m0[49]     9.66
##    m0[50]    47.50
##    y0[1]     -6.47
##    y0[2]     -8.59
##    y0[3]     38.94
##    y0[4]     24.68
##    y0[5]     35.80
##    y0[6]     45.48
##    y0[7]     22.54
##    y0[8]      9.63
##    y0[9]     51.73
##    y0[10]    44.26
##    y0[11]   -22.82
##    y0[12]    54.11
##    y0[13]    60.51
##    y0[14]    22.77
##    y0[15]    33.30
##    y0[16]    23.19
##    y0[17]    19.22
##    y0[18]    18.87
##    y0[19]    16.99
##    y0[20]    21.14
##    y0[21]    -4.88
##    y0[22]    32.99
##    y0[23]    25.37
##    y0[24]    16.04
##    y0[25]    23.39
##    y0[26]    11.39
##    y0[27]    45.55
##    y0[28]    16.25
##    y0[29]    30.36
##    y0[30]    41.61
##    y0[31]    21.84
##    y0[32]    27.37
##    y0[33]     8.73
##    y0[34]    40.02
##    y0[35]    51.81
##    y0[36]    38.49
##    y0[37]    41.32
##    y0[38]    40.55
##    y0[39]    26.07
##    y0[40]    49.59
##    y0[41]    47.88
##    y0[42]    47.04
##    y0[43]     2.00
##    y0[44]    19.31
##    y0[45]    28.59
##    y0[46]     7.95
##    y0[47]    45.51
##    y0[48]    28.33
##    y0[49]    27.21
##    y0[50]    39.47
mcmc=goMCMC(mdl,data,wrm=1000,smp=2000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 2 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 3 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 4 Iteration:    1 / 3000 [  0%]  (Warmup) 
## Chain 1 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 1 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 2 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 2 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 3 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 3 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 4 Iteration: 1000 / 3000 [ 33%]  (Warmup) 
## Chain 4 Iteration: 1001 / 3000 [ 33%]  (Sampling) 
## Chain 1 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 2 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 3 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 4 Iteration: 2000 / 3000 [ 66%]  (Sampling) 
## Chain 1 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 1 finished in 1.2 seconds.
## Chain 2 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 3 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 4 Iteration: 3000 / 3000 [100%]  (Sampling) 
## Chain 2 finished in 1.3 seconds.
## Chain 3 finished in 1.3 seconds.
## Chain 4 finished in 1.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 1.3 seconds.
## Total execution time: 1.4 seconds.
seeMCMC(mcmc,'m')
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   -116.83 -116.60 5.04 4.69 -125.45 -109.01 1.00      830      355
##    b00       7.89    7.90 3.80 3.40    1.85   13.90 1.00     4118     2995
##    sb0       9.57    8.88 3.64 3.04    5.10   16.51 1.00     4730     4995
##    b0[1]    16.28   16.44 7.04 6.97    4.60   27.51 1.00     3612     1665
##    b0[2]    10.98   10.84 4.12 4.08    4.42   17.88 1.00     2449     2785
##    b0[3]    19.46   19.57 2.99 2.94   14.49   24.23 1.00     6108     4781
##    b0[4]     0.36    0.24 3.87 3.81   -5.77    6.80 1.00     5684     3885
##    b0[5]     7.24    7.28 3.95 3.96    0.78   13.71 1.00     2590     1211
##    b0[6]     4.81    4.78 2.31 2.28    1.07    8.81 1.00     2195      672
##    b0[7]     7.28    7.27 4.34 4.18    0.22   14.45 1.00     6415     5961
##    b0[8]    -2.17   -2.27 3.71 3.70   -8.04    4.04 1.00     3033     4060
##    b0[9]     6.77    6.98 4.48 4.34   -1.02   13.66 1.00     4506     4735
##    b10       1.82    1.82 0.26 0.22    1.42    2.23 1.00     5230     5290
##    sb1       0.53    0.47 0.32 0.26    0.13    1.10 1.00      749      295
##    b1[1]     2.17    2.10 0.55 0.49    1.40    3.18 1.00     4375     4999
##    b1[2]     1.37    1.38 0.29 0.30    0.89    1.84 1.00     1515      715
##    b1[3]     1.83    1.83 0.24 0.23    1.45    2.22 1.00     6103     4864
##    b1[4]     1.52    1.61 0.59 0.49    0.39    2.33 1.00     5070     4319
##    b1[5]     2.01    2.01 0.27 0.27    1.58    2.45 1.00     4253     5535
##    b1[6]     2.13    2.13 0.23 0.23    1.75    2.50 1.00     2406     1383
##    b1[7]     1.72    1.73 0.31 0.30    1.20    2.22 1.00     6744     5543
##    b1[8]     1.54    1.56 0.31 0.31    1.00    2.01 1.00     3450     5350
##    b1[9]     2.16    2.10 0.49 0.45    1.47    3.04 1.00     4213     4278
##    s         4.49    4.45 0.57 0.55    3.65    5.49 1.00     4858     5462
##    m0[1]    13.22   13.20 1.69 1.65   10.50   16.14 1.00     2866      951
##    m0[2]    15.96   15.94 2.93 2.95   11.21   20.76 1.00     2379     1251
##    m0[3]    34.06   34.00 2.21 2.18   30.50   37.85 1.00     2083      926
##    m0[4]    28.64   28.71 2.04 1.98   25.24   31.85 1.00     6095     6002
##    m0[5]    46.43   46.41 2.27 2.22   42.78   50.17 1.00     6814     5594
##    m0[6]    37.90   37.87 2.75 2.74   33.43   42.63 1.00     1694      776
##    m0[7]    15.47   15.45 1.71 1.68   12.66   18.27 1.00     7277     4665
##    m0[8]    25.94   26.04 2.30 2.25   22.14   29.58 1.00     6106     5587
##    m0[9]    49.79   49.78 1.95 1.95   46.62   52.96 1.00     6203     5561
##    m0[10]   42.18   42.30 4.36 4.27   34.87   49.18 1.00     5719     4167
##    m0[11]    6.31    6.28 2.18 2.16    2.77   10.10 1.00     2260      731
##    m0[12]   46.06   46.08 2.92 2.89   41.27   50.86 1.00     4855     4103
##    m0[13]   46.28   46.26 2.25 2.20   42.66   50.00 1.00     6816     5533
##    m0[14]   19.34   19.33 1.98 1.97   16.07   22.56 1.00     6974     4742
##    m0[15]   24.04   24.02 1.42 1.40   21.73   26.35 1.00     7464     5450
##    m0[16]   27.82   27.82 1.83 1.83   24.88   30.83 1.00     2318     1381
##    m0[17]   18.61   18.62 2.74 2.74   14.07   22.99 1.00     5611     3492
##    m0[18]    8.50    8.49 2.05 2.02    5.15   11.88 1.00     4283     4903
##    m0[19]    9.16    9.14 1.96 1.93    6.00   12.53 1.00     2424      890
##    m0[20]   16.91   16.82 3.03 3.01   12.03   21.98 1.00     3155     4164
##    m0[21]    3.88    3.84 3.22 3.20   -1.39    9.30 1.00     6481     4079
##    m0[22]   21.17   21.19 2.20 2.20   17.52   24.70 1.00     6348     4749
##    m0[23]   24.83   24.90 2.72 2.73   20.24   29.19 1.00     5634     5314
##    m0[24]   14.63   14.64 3.36 3.27    9.18   20.15 1.00     6068     6436
##    m0[25]   34.91   34.93 1.58 1.56   32.28   37.45 1.00     6085     5966
##    m0[26]   17.07   17.06 1.50 1.48   14.65   19.68 1.00     3864     1208
##    m0[27]   47.45   47.45 1.76 1.76   44.56   50.34 1.00     6186     4117
##    m0[28]   22.58   22.59 1.41 1.38   20.29   24.88 1.00     6814     6080
##    m0[29]   33.12   33.04 3.42 3.33   27.57   38.89 1.00     5424     3682
##    m0[30]   38.41   38.40 1.70 1.67   35.68   41.19 1.00     5188     4868
##    m0[31]   20.32   20.36 2.62 2.62   16.00   24.54 1.00     6011     3603
##    m0[32]   34.52   34.54 1.60 1.57   31.86   37.08 1.00     6094     6039
##    m0[33]    5.93    5.92 2.38 2.37    2.04    9.90 1.00     3740     5018
##    m0[34]   54.08   54.07 2.35 2.33   50.28   57.91 1.00     6288     5774
##    m0[35]   31.56   31.54 1.96 1.93   28.36   34.80 1.00     2776     1263
##    m0[36]   10.02   10.01 3.38 3.30    4.45   15.64 1.00     6434     5095
##    m0[37]   20.00   19.99 2.53 2.50   15.85   24.22 1.00     4018     5127
##    m0[38]   38.82   38.80 1.72 1.68   36.06   41.64 1.00     5346     4916
##    m0[39]   33.68   33.65 1.62 1.61   31.09   36.38 1.00     3026     1346
##    m0[40]   32.27   32.28 2.75 2.68   27.80   36.72 1.00     6690     4214
##    m0[41]   32.72   32.69 1.63 1.61   30.10   35.44 1.00     2837     1350
##    m0[42]   52.02   51.99 2.15 2.15   48.56   55.52 1.00     6250     5641
##    m0[43]   12.43   12.42 1.74 1.70    9.63   15.45 1.00     2748      953
##    m0[44]   27.25   27.26 1.51 1.48   24.78   29.72 1.00     8183     5133
##    m0[45]   34.93   34.88 2.32 2.28   31.20   38.88 1.00     1949      893
##    m0[46]   38.17   38.22 3.28 3.21   32.79   43.49 1.00     6706     4179
##    m0[47]   36.64   36.68 1.51 1.49   34.13   39.05 1.00     6145     6094
##    m0[48]   32.16   32.20 1.75 1.71   29.21   34.93 1.00     6068     6292
##    m0[49]    7.52    7.51 2.17 2.14    3.96   11.10 1.00     4055     4767
##    m0[50]   42.76   42.77 2.62 2.60   38.47   47.04 1.00     5359     4143
##    y0[1]    13.28   13.20 4.81 4.67    5.53   21.28 1.00     8258     8144
##    y0[2]    15.90   15.90 5.35 5.22    6.97   24.77 1.00     6654     7256
##    y0[3]    34.05   34.06 4.99 4.98   25.97   42.21 1.00     7190     7706
##    y0[4]    28.60   28.58 4.93 4.79   20.68   36.75 1.00     7469     7821
##    y0[5]    46.49   46.52 5.05 4.99   38.15   54.72 1.00     7663     7126
##    y0[6]    37.96   37.91 5.22 5.16   29.51   46.68 1.00     6884     6871
##    y0[7]    15.50   15.44 4.86 4.80    7.56   23.44 1.00     8019     7363
##    y0[8]    26.01   26.02 5.00 4.91   17.86   34.15 1.00     7951     7672
##    y0[9]    49.84   49.85 4.89 4.80   41.82   57.77 1.00     7868     7973
##    y0[10]   42.21   42.17 6.29 6.19   32.03   52.67 1.00     8362     7510
##    y0[11]    6.31    6.23 5.01 4.87   -1.83   14.60 1.00     7173     7346
##    y0[12]   46.16   46.16 5.34 5.29   37.42   55.04 1.00     6264     7520
##    y0[13]   46.34   46.39 5.09 5.06   38.05   54.77 1.00     7682     7210
##    y0[14]   19.29   19.27 4.86 4.77   11.36   27.25 1.00     6412     7675
##    y0[15]   23.93   23.94 4.79 4.67   16.07   31.80 1.00     8057     7823
##    y0[16]   27.76   27.79 4.95 4.90   19.70   35.96 1.00     6546     6999
##    y0[17]   18.60   18.59 5.35 5.22    9.82   27.38 1.00     7258     7555
##    y0[18]    8.45    8.48 4.94 4.86    0.17   16.63 1.00     7126     7178
##    y0[19]    9.12    9.12 4.95 4.91    0.94   17.26 1.00     7264     7303
##    y0[20]   16.95   16.99 5.46 5.35    7.92   25.85 1.00     5686     6882
##    y0[21]    3.88    3.83 5.59 5.47   -5.29   12.99 1.00     7259     7484
##    y0[22]   21.20   21.27 4.95 4.88   13.11   29.30 1.00     7226     7259
##    y0[23]   24.78   24.81 5.27 5.14   16.08   33.45 1.00     6567     7769
##    y0[24]   14.64   14.58 5.65 5.60    5.52   24.22 1.00     7575     7449
##    y0[25]   34.94   34.95 4.90 4.83   26.84   42.95 1.00     7661     7374
##    y0[26]   17.12   17.10 4.87 4.72    9.14   24.88 1.00     7516     7726
##    y0[27]   47.50   47.47 4.83 4.74   39.51   55.40 1.00     7561     7589
##    y0[28]   22.57   22.63 4.75 4.69   14.79   30.22 1.00     7934     7486
##    y0[29]   33.06   33.06 5.73 5.65   23.71   42.46 1.00     7563     7209
##    y0[30]   38.42   38.39 4.83 4.73   30.45   46.38 1.00     7657     7721
##    y0[31]   20.35   20.32 5.26 5.14   11.78   29.00 1.00     8028     7635
##    y0[32]   34.55   34.57 4.79 4.61   26.64   42.42 1.00     7674     7872
##    y0[33]    5.97    5.92 5.22 5.11   -2.65   14.57 1.00     7727     7931
##    y0[34]   54.18   54.14 5.10 4.92   45.70   62.58 1.00     7647     7877
##    y0[35]   31.56   31.57 4.89 4.69   23.46   39.64 1.00     6710     7857
##    y0[36]   10.06   10.03 5.67 5.69    0.69   19.38 1.00     8136     6913
##    y0[37]   19.99   19.99 5.18 5.15   11.39   28.35 1.00     7407     7749
##    y0[38]   38.80   38.77 4.78 4.61   30.87   46.63 1.00     7407     7516
##    y0[39]   33.67   33.68 4.80 4.69   25.75   41.66 1.00     7398     7442
##    y0[40]   32.19   32.24 5.28 5.23   23.56   40.85 1.00     8265     7337
##    y0[41]   32.71   32.66 4.81 4.75   24.95   40.65 1.00     7946     7759
##    y0[42]   51.97   51.95 4.94 4.84   43.82   60.08 1.00     8211     8027
##    y0[43]   12.47   12.45 4.82 4.71    4.64   20.34 1.00     7399     7045
##    y0[44]   27.41   27.45 4.76 4.69   19.69   35.05 1.00     8157     7715
##    y0[45]   34.97   35.00 5.06 4.96   26.68   43.26 1.00     6921     7626
##    y0[46]   38.20   38.23 5.65 5.55   28.85   47.50 1.00     7330     6932
##    y0[47]   36.74   36.71 4.73 4.56   28.84   44.56 1.00     7638     7483
##    y0[48]   32.11   32.14 4.85 4.70   24.17   40.09 1.00     7992     7425
##    y0[49]    7.51    7.50 5.05 4.86   -0.82   15.67 1.00     7264     7583
##    y0[50]   42.72   42.71 5.20 5.13   34.21   51.16 1.00     7826     7712

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,c=c,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,shape=c,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  #geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  #geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  #geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

qplot(x,y,facets=~c,size=I(2),col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

ex13

generalized linear mixed model

(X,y)i=1-n
b[b0,b1,...]
linear model    y~N(Xb,s)  
generalized linear model    y~dist.(m=link(Xb),s)  

fixed effect    b0, b1,...  
individual random effect   b0+r0i r0~N(0,sr0), b1+r1i r1~N(0,sr1),...  
class c=1-k  
class effect    b0+r0c r0~N(0,sr0), b1+r1c r1~N(0,sr1),...  

for y=dist.(m,s)
random intercept model    m=b0+r0i+b1*x, m=b0+r0c+b1*x  
random coefficient model  m=b0+(b1+r1i)*x, m=b0+(b1+r1c)*x  
mixed model   m=b0+r0i+(b1+r1i)*x, m=b0+r0c+(b1+r1c)*x  

note  
@ yi=b0+b1*xi+r0i is not useful, ri is included in s  
@ yi=b0+b1*xi+r0c is useful, repeated observation can be treated by class effect
@ when appling Poisson dist.(y is larger than 0, error is larger at large x),
  but variance is larger than mean (over dispersion),
  random intercept model is useful
n=20
x=runif(n,0,10)
r0=rnorm(n,0,1)
y=rpois(n,exp(2+0.1*x+r0))
qplot(x,y)

ex13-0.stan

data{
  int N;
  vector[N] x;
  array[N] int y;
}
parameters{
  real b0;
  real b1;
}
model{
  y~poisson_log(b0+b1*x);
}
generated quantities{
  vector[N] m0;
  array[N] int y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=poisson_log_rng(m0[i]);
  }
}
mdl=cmdstan_model('./ex13-0.stan') 
data=list(N=n,x=x,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = 274.414 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       11       318.903   0.000798363     0.0165867      0.9917      0.9917       14    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     318.90
##    b0         1.90
##    b1         0.10
##    m0[1]      2.68
##    m0[2]      1.91
##    m0[3]      2.31
##    m0[4]      2.12
##    m0[5]      2.71
##    m0[6]      2.05
##    m0[7]      2.61
##    m0[8]      2.81
##    m0[9]      2.81
##    m0[10]     2.12
##    m0[11]     2.18
##    m0[12]     2.50
##    m0[13]     2.45
##    m0[14]     2.05
##    m0[15]     2.04
##    m0[16]     2.77
##    m0[17]     2.33
##    m0[18]     2.37
##    m0[19]     2.29
##    m0[20]     2.14
##    y0[1]     13.00
##    y0[2]      8.00
##    y0[3]     17.00
##    y0[4]      5.00
##    y0[5]      9.00
##    y0[6]      3.00
##    y0[7]     12.00
##    y0[8]     18.00
##    y0[9]     22.00
##    y0[10]    12.00
##    y0[11]     9.00
##    y0[12]    13.00
##    y0[13]    15.00
##    y0[14]    12.00
##    y0[15]     9.00
##    y0[16]    14.00
##    y0[17]     9.00
##    y0[18]     9.00
##    y0[19]     9.00
##    y0[20]    11.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m')
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   318.01 318.31 0.94 0.63 316.16 318.86 1.00      531      603
##    b0       1.90   1.89 0.13 0.13   1.69   2.11 1.01      458      367
##    b1       0.10   0.10 0.02 0.02   0.06   0.13 1.01      568      515
##    m0[1]    2.68   2.68 0.09 0.08   2.53   2.82 1.00     1678     1264
##    m0[2]    1.91   1.91 0.13 0.13   1.70   2.12 1.01      458      372
##    m0[3]    2.31   2.31 0.07 0.07   2.19   2.42 1.01      631      831
##    m0[4]    2.12   2.12 0.09 0.09   1.97   2.28 1.01      473      387
##    m0[5]    2.71   2.71 0.09 0.09   2.56   2.85 1.00     1561     1264
##    m0[6]    2.06   2.05 0.10 0.10   1.89   2.23 1.01      461      375
##    m0[7]    2.60   2.60 0.08 0.07   2.48   2.72 1.00     1970     1301
##    m0[8]    2.80   2.81 0.11 0.11   2.63   2.98 1.00     1254     1313
##    m0[9]    2.80   2.80 0.11 0.11   2.62   2.97 1.00     1269     1296
##    m0[10]   2.12   2.12 0.09 0.09   1.97   2.28 1.01      473      384
##    m0[11]   2.18   2.19 0.08 0.08   2.05   2.33 1.01      497      499
##    m0[12]   2.49   2.50 0.07 0.06   2.38   2.60 1.00     1556     1384
##    m0[13]   2.44   2.44 0.06 0.06   2.34   2.55 1.01     1163     1233
##    m0[14]   2.05   2.05 0.11 0.11   1.88   2.22 1.01      460      375
##    m0[15]   2.04   2.04 0.11 0.11   1.87   2.22 1.01      460      375
##    m0[16]   2.76   2.76 0.10 0.10   2.60   2.92 1.00     1363     1367
##    m0[17]   2.33   2.33 0.07 0.06   2.22   2.44 1.01      677      929
##    m0[18]   2.36   2.36 0.07 0.06   2.25   2.47 1.01      778      978
##    m0[19]   2.29   2.29 0.07 0.07   2.18   2.41 1.01      605      804
##    m0[20]   2.14   2.14 0.09 0.09   1.99   2.29 1.01      478      426
##    y0[1]   14.55  14.00 3.99 4.45   8.00  21.00 1.00     2075     1926
##    y0[2]    6.80   7.00 2.74 2.97   3.00  12.00 1.00     1637     1760
##    y0[3]   10.01  10.00 3.18 2.97   5.00  15.00 1.00     1803     1814
##    y0[4]    8.32   8.00 3.02 2.97   4.00  13.05 1.00     1583     1876
##    y0[5]   14.92  15.00 4.06 4.45   9.00  22.00 1.00     1872     1915
##    y0[6]    7.87   8.00 2.94 2.97   3.00  13.00 1.00     1536     1868
##    y0[7]   13.41  13.00 3.80 2.97   8.00  20.00 1.00     2108     1852
##    y0[8]   16.56  16.00 4.42 4.45  10.00  24.00 1.00     2048     1983
##    y0[9]   16.45  16.00 4.45 4.45  10.00  24.00 1.00     1930     1988
##    y0[10]   8.41   8.00 2.99 2.97   4.00  14.00 1.00     1885     1877
##    y0[11]   9.06   9.00 3.08 2.97   4.00  14.00 1.00     1721     1751
##    y0[12]  12.16  12.00 3.68 2.97   6.00  18.00 1.00     2089     1828
##    y0[13]  11.46  11.00 3.39 2.97   6.00  17.00 1.00     2073     2037
##    y0[14]   7.78   8.00 2.87 2.97   3.00  13.00 1.00     1746     2097
##    y0[15]   7.74   8.00 2.89 2.97   3.00  13.00 1.00     1683     1964
##    y0[16]  16.04  16.00 4.36 4.45   9.00  24.00 1.00     1944     1867
##    y0[17]  10.21  10.00 3.35 2.97   5.00  16.00 1.00     1804     1870
##    y0[18]  10.79  11.00 3.44 2.97   5.00  17.00 1.00     1887     1949
##    y0[19]   9.87  10.00 3.10 2.97   5.00  15.00 1.00     2014     1997
##    y0[20]   8.51   8.00 3.04 2.97   4.00  14.00 1.00     1736     1890

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=exp(m)),xy,col='black')

ex13-1.stan

data{
  int N;
  vector[N] x;
  array[N] int y;
}
parameters{
  real b0;
  real b1;
  real<lower=0> sr0;
  vector[N] r0;
}
model{
  r0~normal(0,sr0);
  for(i in 1:N)
    y~poisson_log(b0+r0[i]+b1*x);
}
generated quantities{
  vector[N] m0;
  array[N] int y0;
  for(i in 1:N){
    m0[i]=b0+r0[i]+b1*x[i];
    y0[i]=poisson_log_rng(m0[i]);
  }
}
mdl=cmdstan_model('./ex13-1.stan') 
data=list(N=n,x=x,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -3847.38 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       99       6539.77     0.0415256       8906.86      0.2214      0.2214      146    
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##      152       6658.97    0.00244894   1.88547e+06      0.8696      0.8696      246    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__    6658.97
##    b0         1.85
##    b1         0.14
##    sr0        0.00
##    r0[1]      0.00
##    r0[2]      0.00
##    r0[3]      0.00
##    r0[4]      0.00
##    r0[5]      0.00
##    r0[6]      0.00
##    r0[7]      0.00
##    r0[8]      0.00
##    r0[9]      0.00
##    r0[10]     0.00
##    r0[11]     0.00
##    r0[12]     0.00
##    r0[13]     0.00
##    r0[14]     0.00
##    r0[15]     0.00
##    r0[16]     0.00
##    r0[17]     0.00
##    r0[18]     0.00
##    r0[19]     0.00
##    r0[20]     0.00
##    m0[1]      2.98
##    m0[2]      1.87
##    m0[3]      2.44
##    m0[4]      2.17
##    m0[5]      3.03
##    m0[6]      2.08
##    m0[7]      2.87
##    m0[8]      3.17
##    m0[9]      3.16
##    m0[10]     2.18
##    m0[11]     2.26
##    m0[12]     2.72
##    m0[13]     2.64
##    m0[14]     2.07
##    m0[15]     2.05
##    m0[16]     3.11
##    m0[17]     2.47
##    m0[18]     2.53
##    m0[19]     2.42
##    m0[20]     2.20
##    y0[1]     24.00
##    y0[2]      4.00
##    y0[3]     13.00
##    y0[4]     10.00
##    y0[5]     25.00
##    y0[6]      5.00
##    y0[7]     16.00
##    y0[8]     22.00
##    y0[9]     20.00
##    y0[10]    12.00
##    y0[11]     2.00
##    y0[12]    14.00
##    y0[13]    15.00
##    y0[14]     6.00
##    y0[15]     6.00
##    y0[16]    20.00
##    y0[17]    12.00
##    y0[18]    14.00
##    y0[19]    12.00
##    y0[20]    15.00
mcmc=goMCMC(mdl,data,wrm=500)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
## Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
## Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 2 finished in 0.7 seconds.
## Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 3 finished in 0.8 seconds.
## Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
## Chain 1 finished in 0.9 seconds.
## Chain 4 finished in 0.8 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.8 seconds.
## Total execution time: 1.0 seconds.
seeMCMC(mcmc,'[m,r]')
##  variable    mean  median    sd   mad      q5     q95 rhat ess_bulk ess_tail
##    lp__   6454.08 6452.80 14.62 17.20 6431.54 6477.87 1.06       61       96
##    b0        1.90    1.90  0.03  0.03    1.85    1.95 1.01      444      913
##    b1        0.10    0.10  0.01  0.01    0.09    0.11 1.01      424     1034
##    sr0       0.01    0.01  0.01  0.01    0.00    0.03 1.05       63      111
##    r0[1]     0.00    0.00  0.01  0.01   -0.02    0.03 1.02     2753      857
##    r0[2]     0.00    0.00  0.02  0.01   -0.03    0.02 1.01     3486      863
##    r0[3]     0.00    0.00  0.02  0.01   -0.02    0.03 1.02     3180      802
##    r0[4]     0.00    0.00  0.02  0.01   -0.02    0.02 1.02     4061      452
##    r0[5]     0.00    0.00  0.01  0.01   -0.02    0.02 1.02     3127      838
##    r0[6]     0.00    0.00  0.02  0.01   -0.03    0.02 1.02     2805      720
##    r0[7]     0.00    0.00  0.02  0.01   -0.03    0.02 1.04     3817      799
##    r0[8]     0.00    0.00  0.01  0.01   -0.02    0.02 1.03     3862      646
##    r0[9]     0.00    0.00  0.02  0.01   -0.03    0.03 1.04     2789      860
##    r0[10]    0.00    0.00  0.02  0.01   -0.03    0.03 1.02     4005      832
##    r0[11]    0.00    0.00  0.02  0.01   -0.02    0.02 1.03     3282      733
##    r0[12]    0.00    0.00  0.02  0.01   -0.03    0.03 1.02     3616      784
##    r0[13]    0.00    0.00  0.02  0.01   -0.03    0.02 1.02     3440      630
##    r0[14]    0.00    0.00  0.02  0.01   -0.03    0.03 1.03     3249      751
##    r0[15]    0.00    0.00  0.02  0.01   -0.02    0.02 1.02     3095      919
##    r0[16]    0.00    0.00  0.01  0.01   -0.02    0.02 1.01     2387      508
##    r0[17]    0.00    0.00  0.02  0.01   -0.03    0.02 1.03     3014      750
##    r0[18]    0.00    0.00  0.02  0.01   -0.02    0.02 1.02     3729      624
##    r0[19]    0.00    0.00  0.01  0.01   -0.02    0.02 1.02     3096      934
##    r0[20]    0.00    0.00  0.02  0.01   -0.02    0.02 1.03     4177      558
##    m0[1]     2.68    2.68  0.02  0.02    2.64    2.72 1.00     1961     1407
##    m0[2]     1.91    1.91  0.04  0.04    1.85    1.96 1.01      519      868
##    m0[3]     2.31    2.31  0.02  0.02    2.27    2.34 1.00     1491     1371
##    m0[4]     2.12    2.12  0.03  0.03    2.08    2.17 1.01      735      852
##    m0[5]     2.71    2.71  0.02  0.02    2.67    2.75 1.00     1292     1573
##    m0[6]     2.05    2.05  0.03  0.03    2.01    2.10 1.00      613     1058
##    m0[7]     2.61    2.61  0.02  0.02    2.57    2.64 1.02     2856     1257
##    m0[8]     2.81    2.81  0.03  0.03    2.77    2.86 1.00     1072     1543
##    m0[9]     2.81    2.81  0.03  0.03    2.76    2.85 1.00     1062     1568
##    m0[10]    2.12    2.12  0.03  0.03    2.08    2.17 1.00      811     1492
##    m0[11]    2.18    2.18  0.03  0.03    2.14    2.22 1.00      930     1025
##    m0[12]    2.50    2.50  0.02  0.02    2.46    2.53 1.01     2913     1250
##    m0[13]    2.45    2.45  0.02  0.02    2.41    2.48 1.02     2630     1082
##    m0[14]    2.05    2.05  0.03  0.03    2.00    2.10 1.01      734     1025
##    m0[15]    2.04    2.04  0.03  0.03    1.99    2.09 1.00      555     1289
##    m0[16]    2.77    2.77  0.03  0.03    2.73    2.81 1.00     1215     1528
##    m0[17]    2.33    2.33  0.02  0.02    2.29    2.36 1.01     1604     1008
##    m0[18]    2.37    2.37  0.02  0.02    2.33    2.40 1.02     1982     1217
##    m0[19]    2.29    2.29  0.02  0.02    2.26    2.33 1.00     1130     1066
##    m0[20]    2.14    2.14  0.03  0.03    2.09    2.18 1.00      816     1343
##    y0[1]    14.71   15.00  3.81  4.45    9.00   21.00 1.00     1711     1469
##    y0[2]     6.77    7.00  2.57  2.97    3.00   11.00 1.00     1734     1827
##    y0[3]     9.99   10.00  3.16  2.97    5.00   16.00 1.00     2016     1872
##    y0[4]     8.28    8.00  2.78  2.97    4.00   13.00 1.00     2116     1846
##    y0[5]    15.10   15.00  3.95  4.45    9.00   22.00 1.00     1915     1944
##    y0[6]     7.93    8.00  2.83  2.97    3.95   13.00 1.00     1944     1911
##    y0[7]    13.46   13.00  3.58  2.97    8.00   20.00 1.00     2058     2014
##    y0[8]    16.65   17.00  4.06  4.45   10.00   23.05 1.00     1976     2053
##    y0[9]    16.54   16.00  4.09  4.45   10.00   23.00 1.00     1951     1958
##    y0[10]    8.36    8.00  2.88  2.97    4.00   13.00 1.00     2056     1983
##    y0[11]    8.83    9.00  2.97  2.97    4.00   14.00 1.00     2017     1821
##    y0[12]   12.17   12.00  3.58  2.97    6.00   18.00 1.00     1853     2006
##    y0[13]   11.50   11.00  3.45  2.97    6.00   18.00 1.00     1900     1983
##    y0[14]    7.72    8.00  2.76  2.97    4.00   13.00 1.00     2114     1902
##    y0[15]    7.82    8.00  2.84  2.97    3.00   13.00 1.00     2089     1977
##    y0[16]   15.89   16.00  3.93  4.45   10.00   23.00 1.00     2060     1680
##    y0[17]   10.32   10.00  3.24  2.97    5.00   16.00 1.00     1936     2068
##    y0[18]   10.64   10.50  3.26  3.71    6.00   16.00 1.00     1922     1898
##    y0[19]    9.99   10.00  3.19  2.97    5.00   15.00 1.00     1935     1899
##    y0[20]    8.43    8.00  2.85  2.97    4.00   13.00 1.00     1866     1805

m0=mcmc$draws('m0')
smm0=summary(m0)

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

xy=tibble(x=x,m=smm0$median,ml5=smm0$q5,mu5=smm0$q95,yl5=smy0$q5,yu5=smy0$q95)

qplot(x,y,size=I(2),col=I('red'))+
  scale_shape_manual(values=1:k) + 
  geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=exp(m)),xy,col='black')